##########################################################
# John Henderson and Alex Theodoridis
# Replication Data for: "Seeing Spots", 
#  Forthcoming in Political Behavior, August 20, 2017
# 
##########################################################
#
#  table_ii_vii.R
#  -- file reproduces regression results for tables ii - vii in the appendix
#
##########################################################
    
rm(list=ls())  
set.seed(1005)                                  
source('~/Dropbox/Seeing_Spots/replication/covar_data.R')  # contains pre_data

ymat=cbind(video_skipped, replay, share, getlink, all_y,time_watched, total_time)
                        
# SUR for coherency interpretations of coefficients across all outcomes (alt. additive outcome)
library(systemfit)
ymat[,1]=1-ymat[,1]

########################################################### 
# experimental regression results + observational controls
########################################################## 
  
party_tr=as.numeric(obama_ad)
party_tr[which(party_tr==1)]=-1
party_tr[which(party_tr==0)]=1

tone_tr=as.numeric(pos_ad)
tone_tr[which(tone_tr==1)]=1
tone_tr[which(tone_tr==0)]=-1

# some interactions and additional variables for below               
news_tone=newsint*tone_tr
women_tone=women*tone_tr 
                                             
# already seen ad instrument
seen_ad=as.numeric(cces_data$AGT308=='Yes')
seen_ad[which(as.numeric(cces_data$AGT308)>3)]=NA
seen_ad_mean=tapply(seen_ad,as.character(tr),mean,na.rm=T)

seen_ad_indx=array(NA,length(seen_ad))
for(i in 1:length(seen_ad_mean)){ 
	ix=which(names(seen_ad_mean)[i]==as.character(tr))	
	seen_ad_indx[ix]=seen_ad_mean[i]
}

# ad is interesting instrument
ad_interesting=as.numeric(cces_data$AGT315_3)
ad_interesting[which(ad_interesting==4)]=NA
ad_interesting=-ad_interesting+2
ad_interesting_mean=tapply(ad_interesting,as.character(tr),mean,na.rm=T)

ad_interesting_indx=array(NA,length(ad_interesting))
for(i in 1:length(ad_interesting_mean)){ 
	ix=which(names(ad_interesting_mean)[i]==as.character(tr))	
	ad_interesting_indx[ix]=ad_interesting_mean[i]
}


###################
# TABLE II - baseline ols

r1 <- ymat[,1]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone
	
r2 <- ymat[,2]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone

r3 <- ymat[,3]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone

r4 <- ymat[,4]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone

fitsur <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='SUR')
summary(fitsur)  

fitols <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='OLS')
summary(fitols)

library(car)
# joint F over SUR
lst=rownames(summary(lm(r4))$coef)   
joint_p=as.data.frame(matrix(1,length(lst),1))

for(j in 1:nrow(joint_p)){ 
	restriction <- c(
		paste("skipreg_",lst[j],"=0",sep=''),
		paste("replayreg_",lst[j],"=0",sep=''),
		paste("sharereg_",lst[j],"=0",sep=''),
		paste("linkreg_",lst[j],"=0",sep=''))
	x=linearHypothesis(fitols, hypothesis.matrix=restriction, test = "F")
	joint_p[j,1]=(x[[4]])[2] 
	rownames(joint_p)[j]=lst[j]
}
 
t(t(summary(fitsur)$coefficients[,1]))
t(t(summary(fitsur)$coefficients[,2]))
t(t(summary(fitols)$coefficients[,1]))
t(t(summary(fitols)$coefficients[,2]))
round(joint_p,digits=3)
           
###################
# TABLE III - logit     

log1=summary(glm(r1,family=binomial(link=logit)))
log2=summary(glm(r2,family=binomial(link=logit)))
log3=summary(glm(r3,family=binomial(link=logit)))
log4=summary(glm(r4,family=binomial(link=logit))) 
  
###################
# TABLE IV - interaction of tone x interesting ad
  

r1 <- ymat[,1]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+ad_interesting_indx*tone_tr
   

r2 <- ymat[,2]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+ad_interesting_indx*tone_tr


r3 <- ymat[,3]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+ad_interesting_indx*tone_tr


r4 <- ymat[,4]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+ad_interesting_indx*tone_tr
  
fitsur <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='SUR')
summary(fitsur)

fitols <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='OLS')
summary(fitols)               
   
# skipping                                                             
 #ad_interesting_indx          0.10901723  0.40700049  0.26786 0.7888643 
 #tone_tr:ad_interesting_indx  0.54696266  0.40631863  1.34614 0.1785526  
  
library(car)
# joint F over SUR
lst=rownames(summary(lm(r4))$coef)   
joint_p=as.data.frame(matrix(1,length(lst),1))

for(j in 1:nrow(joint_p)){ 
	restriction <- c(
		paste("skipreg_",lst[j],"=0",sep=''),
		paste("replayreg_",lst[j],"=0",sep=''),
		paste("sharereg_",lst[j],"=0",sep=''),
		paste("linkreg_",lst[j],"=0",sep=''))
	x=linearHypothesis(fitols, hypothesis.matrix=restriction, test = "F")
	joint_p[j,1]=(x[[4]])[2] 
	rownames(joint_p)[j]=lst[j]
}

t(t(summary(fitsur)$coefficients[,1]))
t(t(summary(fitsur)$coefficients[,2]))
t(t(summary(fitols)$coefficients[,1]))
t(t(summary(fitols)$coefficients[,2]))


###################
# TABLE V - interaction of tone x seen ad

r1 <- ymat[,1]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+seen_ad_indx*tone_tr
   

r2 <- ymat[,2]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+seen_ad_indx*tone_tr


r3 <- ymat[,3]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+seen_ad_indx*tone_tr


r4 <- ymat[,4]~
party_tr*pid_lean*tone_tr+
newsint+pid_strenth+  
reg+uncertain+battleground+
abs(mdsscale)+
women+black+educ+income+married+
unemployed+
vote+churchattend+news_tone+seen_ad_indx*tone_tr
   
fitsur <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='SUR')
summary(fitsur)

fitols <- systemfit(list(skipreg = r1, replayreg = r2, sharereg = r3, linkreg = r4),method='OLS')
summary(fitols)               
   
# skipping                                                             
 #seen_ad_indx              -0.31944681  0.18224018 -1.75289  0.0799183 .  
 #tone_tr:seen_ad_indx       0.18758694  0.18211677  1.03004  0.3032344    

library(car)
# joint F over SUR
lst=rownames(summary(lm(r4))$coef)   
joint_p=as.data.frame(matrix(1,length(lst),1))

for(j in 1:nrow(joint_p)){ 
	restriction <- c(
		paste("skipreg_",lst[j],"=0",sep=''),
		paste("replayreg_",lst[j],"=0",sep=''),
		paste("sharereg_",lst[j],"=0",sep=''),
		paste("linkreg_",lst[j],"=0",sep=''))
	x=linearHypothesis(fitols, hypothesis.matrix=restriction, test = "F")
	joint_p[j,1]=(x[[4]])[2] 
	rownames(joint_p)[j]=lst[j]
}

t(t(summary(fitsur)$coefficients[,1]))
t(t(summary(fitsur)$coefficients[,2]))
t(t(summary(fitols)$coefficients[,1]))
t(t(summary(fitols)$coefficients[,2]))


###################
# TABLE VI - post-treatment on skippers

# results for those skipping and not skipping .... 
y1=ymat[,1]
y2=ymat[,2]
y3=ymat[,3]
y4=ymat[,4]

  
r2 <- y2[which(video_skipped==1)]~
party_tr[which(video_skipped==1)]*pid_lean[which(video_skipped==1)]*tone_tr[which(video_skipped==1)]+
newsint[which(video_skipped==1)]+pid_strenth[which(video_skipped==1)]+  
reg[which(video_skipped==1)]+uncertain[which(video_skipped==1)]+battleground[which(video_skipped==1)]+
abs(mdsscale)[which(video_skipped==1)]+#approval+
women[which(video_skipped==1)]+black[which(video_skipped==1)]+educ[which(video_skipped==1)]+income[which(video_skipped==1)]+married[which(video_skipped==1)]+
unemployed[which(video_skipped==1)]+#union+
vote[which(video_skipped==1)]+churchattend[which(video_skipped==1)]+news_tone[which(video_skipped==1)]

r3 <- y3[which(video_skipped==1)]~
party_tr[which(video_skipped==1)]*pid_lean[which(video_skipped==1)]*tone_tr[which(video_skipped==1)]+
newsint[which(video_skipped==1)]+pid_strenth[which(video_skipped==1)]+  
reg[which(video_skipped==1)]+uncertain[which(video_skipped==1)]+battleground[which(video_skipped==1)]+
abs(mdsscale)[which(video_skipped==1)]+#approval+
women[which(video_skipped==1)]+black[which(video_skipped==1)]+educ[which(video_skipped==1)]+income[which(video_skipped==1)]+married[which(video_skipped==1)]+
unemployed[which(video_skipped==1)]+#union+
vote[which(video_skipped==1)]+churchattend[which(video_skipped==1)]+news_tone[which(video_skipped==1)]

r4 <- y4[which(video_skipped==1)]~
party_tr[which(video_skipped==1)]*pid_lean[which(video_skipped==1)]*tone_tr[which(video_skipped==1)]+
newsint[which(video_skipped==1)]+pid_strenth[which(video_skipped==1)]+  
reg[which(video_skipped==1)]+uncertain[which(video_skipped==1)]+battleground[which(video_skipped==1)]+
abs(mdsscale)[which(video_skipped==1)]+#approval+
women[which(video_skipped==1)]+black[which(video_skipped==1)]+educ[which(video_skipped==1)]+income[which(video_skipped==1)]+married[which(video_skipped==1)]+
unemployed[which(video_skipped==1)]+#union+
vote[which(video_skipped==1)]+churchattend[which(video_skipped==1)]+news_tone[which(video_skipped==1)]


post2a=summary(lm(r2))
#459
post3a=summary(lm(r3))
#459
post4a=summary(lm(r4))
#459
     
fitols_1 <- systemfit(list(replayreg = r2, sharereg = r3, linkreg = r4),method='OLS')
summary(fitols_1)

library(car)
# joint F over SUR
lst=rownames(summary(lm(r4))$coef)   
joint_p=as.data.frame(matrix(1,length(lst),1))

for(j in 1:nrow(joint_p)){ 
	restriction <- c(
		#paste("skipreg_",lst[j],"=0",sep=''),
		paste("replayreg_",lst[j],"=0",sep=''),
		paste("sharereg_",lst[j],"=0",sep=''),
		paste("linkreg_",lst[j],"=0",sep=''))
	x=linearHypothesis(fitols_1, hypothesis.matrix=restriction, test = "F")
	joint_p[j,1]=(x[[4]])[2] 
	rownames(joint_p)[j]=lst[j]
}
 
###################
# TABLE VII - post-treatment on non-skippers    
   
r2 <- y2[which(video_skipped==0)]~
party_tr[which(video_skipped==0)]*pid_lean[which(video_skipped==0)]*tone_tr[which(video_skipped==0)]+
newsint[which(video_skipped==0)]+pid_strenth[which(video_skipped==0)]+  
reg[which(video_skipped==0)]+uncertain[which(video_skipped==0)]+battleground[which(video_skipped==0)]+
abs(mdsscale)[which(video_skipped==0)]+#approval+
women[which(video_skipped==0)]+black[which(video_skipped==0)]+educ[which(video_skipped==0)]+income[which(video_skipped==0)]+married[which(video_skipped==0)]+
unemployed[which(video_skipped==0)]+#union+
vote[which(video_skipped==0)]+churchattend[which(video_skipped==0)]+news_tone[which(video_skipped==0)]

r3 <- y3[which(video_skipped==0)]~
party_tr[which(video_skipped==0)]*pid_lean[which(video_skipped==0)]*tone_tr[which(video_skipped==0)]+
newsint[which(video_skipped==0)]+pid_strenth[which(video_skipped==0)]+  
reg[which(video_skipped==0)]+uncertain[which(video_skipped==0)]+battleground[which(video_skipped==0)]+
abs(mdsscale)[which(video_skipped==0)]+#approval+
women[which(video_skipped==0)]+black[which(video_skipped==0)]+educ[which(video_skipped==0)]+income[which(video_skipped==0)]+married[which(video_skipped==0)]+
unemployed[which(video_skipped==0)]+#union+
vote[which(video_skipped==0)]+churchattend[which(video_skipped==0)]+news_tone[which(video_skipped==0)]

r4 <- y4[which(video_skipped==0)]~
party_tr[which(video_skipped==0)]*pid_lean[which(video_skipped==0)]*tone_tr[which(video_skipped==0)]+
newsint[which(video_skipped==0)]+pid_strenth[which(video_skipped==0)]+  
reg[which(video_skipped==0)]+uncertain[which(video_skipped==0)]+battleground[which(video_skipped==0)]+
abs(mdsscale)[which(video_skipped==0)]+#approval+
women[which(video_skipped==0)]+black[which(video_skipped==0)]+educ[which(video_skipped==0)]+income[which(video_skipped==0)]+married[which(video_skipped==0)]+
unemployed[which(video_skipped==0)]+#union+
vote[which(video_skipped==0)]+churchattend[which(video_skipped==0)]+news_tone[which(video_skipped==0)]

 
post2b=summary(lm(r2))
#592
post3b=summary(lm(r3))
#593 
post4b=summary(lm(r4))
#594


fitols_0<- systemfit(list(replayreg = r2, sharereg = r3, linkreg = r4),method='OLS')
summary(fitols_0)

library(car)
# joint F over SUR
lst=rownames(summary(lm(r4))$coef)   
joint_p=as.data.frame(matrix(1,length(lst),1))

for(j in 1:nrow(joint_p)){ 
	restriction <- c(
		#paste("skipreg_",lst[j],"=0",sep=''),
		paste("replayreg_",lst[j],"=0",sep=''),
		paste("sharereg_",lst[j],"=0",sep=''),
		paste("linkreg_",lst[j],"=0",sep=''))
	x=linearHypothesis(fitols_0, hypothesis.matrix=restriction, test = "F")
	joint_p[j,1]=(x[[4]])[2] 
	rownames(joint_p)[j]=lst[j]
}


#END table_ii_vii.R